/*****************************************************************************
*
* Copyright (c) 2000 - 2018, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

// ************************************************************************* //
//                            avtMDSCHISMFileFormat.C                           //
// ************************************************************************* //

#include <avtMDSCHISMFileFormat.h>

#include <string>
#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <avtDatabaseMetaData.h>
#include <DBOptionsAttributes.h>
#include <Expression.h>
#include <InvalidVariableException.h>
#include "FileFormatFavorFactory.h"

#include <vtkCharArray.h>
#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkIntArray.h>
#include <vtkLongArray.h>
#include <vtkShortArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnsignedIntArray.h>
#include <vtkUnsignedShortArray.h>

#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkFieldData.h>
#include <vtkInformation.h>
#include <vtkRectilinearGrid.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <avtGhostData.h>
#include <DebugStream.h>
#ifdef PARALLEL
#include <avtParallel.h>
#endif
using     std::string;
const std::string NODE = MeshConstants10::NODE;
const std::string FACE = MeshConstants10::ELEM;
const std::string SIDE = MeshConstants10::EDGE;

// ****************************************************************************
//  Method: avtMDSCHISMFileFormat constructor
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

avtMDSCHISMFileFormat::avtMDSCHISMFileFormat(const char *filename)
    : avtMTMDFileFormat(filename),
	  m_data_file(filename),
	m_mesh_3d("3D_Mesh"),
	m_layer_mesh("Layer_Mesh"),
	m_mesh_2d("2D_Mesh"),
	m_mesh_2d_no_wet_dry("2D_Mesh_no_wet_dry"),
	m_side_center_point_3d_mesh("side_center_3D"),
	m_side_center_point_2d_mesh("side_center_2D"),
	m_face_center_point_3d_mesh("face_center_3D"),
	m_surface_state_suffix("_surface"),
	m_bottom_state_suffix("_near_bottom"),
	m_depth_average_suffix("_depth_average"),
	m_node_depth_label(MeshConstants10::NODE_DEPTH),
	m_node_depth(MeshConstants10::NODE_DEPTH)
{
    // INITIALIZE DATA MEMBERS
	m_impl = new avtMDSCHISMFileFormatImpl();
	//m_impl.reset(FileFormatFavorFactory::Instance()->CreateInstance("md_nc4"));
	//
	//debug1 << "impl reset to mdnc4 "<<m_impl<<"\n";
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

int
avtMDSCHISMFileFormat::GetNTimesteps(void)
{
	
    return m_impl->GetNTimesteps(m_data_file);
}
void avtMDSCHISMFileFormat::GetTimes(std::vector<double> &a_time)
{
	
	m_impl->GetTimes(a_time);

}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

void
avtMDSCHISMFileFormat::FreeUpResources(void)
{
	 m_impl->FreeUpResources();
}


void
avtMDSCHISMFileFormat::ActivateTimestep(void)
{
	return;
	//m_impl->ActivateTimestep(m_data_file);
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

void
avtMDSCHISMFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{

	
	//
	m_var_name_label_map[m_node_depth_label] = m_node_depth;
	int nblocks = m_impl->num_domain();
	MDSchismOutput * file_ptr = m_impl->get_a_data_file();
	MDSCHISMMeshProvider * mesh_ptr = m_impl->get_a_mesh_provider();

	m_data = file_ptr;
	m_mesh = mesh_ptr;


	avtCentering  nodeCent = AVT_NODECENT;
	avtCentering  zoneCent = AVT_ZONECENT;
	string mesh_name = m_mesh_3d;

	// AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
	// AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
	avtMeshType mt = AVT_UNSTRUCTURED_MESH;

	int block_origin = 0;
	int cell_origin = 0;
	int group_origin = 0;
	int spatial_dimension = 3;
	int topological_dimension = 3;
	double *extents = NULL;
	int * bounds = NULL;
	//add node center 3d mesh
	if (m_mesh->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
	    md->Add(mmd);
	}

	//
   // add layered 2d mesh
	mesh_name = m_layer_mesh;
	topological_dimension = 2;

	if (m_mesh->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		md->Add(mmd);
	}

	//
	// add surface 2d mesh
	//
	mesh_name = m_mesh_2d;

	spatial_dimension = 2;
	topological_dimension = 2;

	avtMeshMetaData *mmd1 = new avtMeshMetaData(bounds, extents, mesh_name,
		nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
	md->Add(mmd1);

	//add 3d side center point mesh
	mesh_name = m_side_center_point_3d_mesh;

	spatial_dimension = 3;
	topological_dimension = 0;

	if (m_mesh->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		md->Add(mmd);
	}

	//add 3d  face center point mesh
	mesh_name = m_face_center_point_3d_mesh;

	spatial_dimension = 3;
	topological_dimension = 0;

	if (m_mesh->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		md->Add(mmd);
	}

	// add 2d side center point mesh
	mesh_name = m_side_center_point_2d_mesh;

	spatial_dimension = 2;
	topological_dimension = 0;

	avtMeshMetaData *mmd2 = new avtMeshMetaData(bounds, extents, mesh_name,
		nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
	md->Add(mmd2);

	//
	// add water surface and depth scalar
	string mesh = m_mesh_2d;
	avtScalarMetaData *smd1 = new avtScalarMetaData(m_node_depth_label, mesh, nodeCent);
	md->Add(smd1);
	// m_var_mesh_map[m_node_surface_label] = mesh;
	m_var_mesh_map[m_node_depth_label] = mesh;

	avtScalarMetaData *smd_node_id = new avtScalarMetaData(MeshConstants10::NODE_GLOBAL_ID, mesh, nodeCent);
	md->Add(smd_node_id);
	m_var_mesh_map[MeshConstants10::NODE_GLOBAL_ID] = mesh;

	avtScalarMetaData *smd_ele_id = new avtScalarMetaData(MeshConstants10::ELE_GLOBAL_ID, mesh, zoneCent);
	md->Add(smd_ele_id);
	m_var_mesh_map[MeshConstants10::ELE_GLOBAL_ID] = mesh;

	//add 3D node level label
	string mesh3d = m_mesh_3d;
	avtScalarMetaData *smd2 = new avtScalarMetaData(MeshConstants10::NODE_LEVEL, mesh3d, nodeCent);
	md->Add(smd2);

	//add 3D element level label
	mesh3d = m_layer_mesh;
	avtScalarMetaData *smd3 = new avtScalarMetaData(MeshConstants10::ELE_LEVEL, mesh3d, zoneCent);
	md->Add(smd3);

	//add 3D side level label
	mesh3d = m_side_center_point_3d_mesh;
	avtScalarMetaData *smd4 = new avtScalarMetaData(MeshConstants10::SIDE_LEVEL, mesh3d, nodeCent);
	md->Add(smd4);

	//add 3D layer lable (for prism center data)
	mesh3d = m_mesh_3d;
	avtScalarMetaData *smd5 = new avtScalarMetaData(MeshConstants10::LAYER, mesh3d, zoneCent);
	md->Add(smd5);

	
	PopulateStateMetaData(md, timeState);
	m_impl->set_var_horizontal_center_map(m_var_horizontal_center_map);
	m_impl->set_var_mesh_map(m_var_mesh_map);
	m_impl->set_var_name_label_map(m_var_name_label_map);
	m_impl->set_var_vertical_center_map(m_var_vertical_center_map);
	return;
	//avtDatabaseMetaData a_md;
    //m_impl->PopulateDatabaseMetaData(&a_md, m_data_file,timeState);
	//*md = a_md;
	//debug1 << "done impl md populate\n";
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::PopulateStateMetaData
//
//  Purpose:
//      Scan the data file and find out all the state variables, 
//      store the name, mesh, and centering information into
//      the input metadata set.
//
//  Programmer: qshu
//  Creation:   Wed Aug 30 08:11:04 PDT 2012
//
// ****************************************************************************

void avtMDSCHISMFileFormat::PopulateStateMetaData(avtDatabaseMetaData * a_metaData,
	int                   a_timeState)
{
	int numVar = m_data->num_vars();
	//debug1 << "get vars " << numVar << endl;
	for (int iVar = 0; iVar < numVar; iVar++)
	{
		//debug1 << iVar;
		SCHISMVar10*  varPtr = m_data->get_var(iVar);
		std::string varName = varPtr->name();

		if (m_data->none_data_var(varName))
		{
			debug1 << varName << " is skipped\n";
			continue;
		}
		debug1 << iVar << " " << varPtr->num_dims() << " " << varName << endl;
		if (!(varPtr->is_defined_over_grid()))
		{
			debug1 << varName << " passed \n";
			continue;
		}

		if (varPtr->is_SCHISM_mesh_parameter())
		{
			debug1 << varName << "passed 2\n";
			continue;
		}

		//debug1 << varName<<"to be added" << endl;

		std::string  location(NODE);
		avtCentering avtCenter(AVT_NODECENT);


		if ((varName == m_node_surface) || (varName == m_node_depth))
		{
			continue;
		}
		std::string  label;
		label = varName;
		// this dic make it easy to find out data set for a visit plot variable
		m_var_name_label_map[label] = varName;

		// handle different for face and node center data
		location = varPtr->get_horizontal_center();

		std::string vertical_center = varPtr->get_vertical_center();

		m_var_horizontal_center_map[varName] = location;
		m_var_vertical_center_map[varName] = vertical_center;

		if (location == FACE)
		{
			//debug1 << " begin add face var meta\n";
			addFaceCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		else if (location == NODE)
		{
			//debug1 << " begin add node var meta\n";
			addNodeCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		else if (location == SIDE)
		{
			//debug1 << " begin add side var meta\n";
			addSideCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		// omit unkown center data
		else
		{
			continue;
		}
	}
}

void    avtMDSCHISMFileFormat::addFaceCenterData(avtDatabaseMetaData * a_metaData,
	SCHISMVar10           * a_varPtr,
	const std::string   & a_varName,
	const std::string   & a_varLabel,
	const avtCentering  & a_center)

{
	// only add face centered  var now 
	string mesh2d = m_mesh_2d;
	string mesh3d = m_mesh_3d;
	std::string level_center = a_varPtr->get_vertical_center();

	if (level_center == MeshConstants10::FULL_LAYER)
	{
		mesh3d = m_layer_mesh;
	}

	avtCentering  faceCent(AVT_ZONECENT);
	// scalar data 2d mesh
	if (a_varPtr->num_dims() <= 2)
	{
		avtScalarMetaData *smd = new avtScalarMetaData(a_varLabel, mesh2d, faceCent);
		a_metaData->Add(smd);

		m_var_mesh_map[a_varLabel] = mesh2d;
		m_var_dim[a_varName] = a_varPtr->num_dims();
	}
	// vector data
	else if (m_impl->SCHISMVarIs3D(a_varPtr))
	{
		if (!(m_mesh->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		if (m_impl->SCHISMVarIsVector(a_varPtr)) // 3d vector
		{
			// last dim is vector component         
			SCHISMDim10* comDim = a_varPtr->get_dim(3);
			int ncomps = comDim->size();
			int ucomps = (ncomps == 2 ? 3 : ncomps);

			avtVectorMetaData *vmd = new avtVectorMetaData(a_varLabel, mesh3d, faceCent, ucomps);
			a_metaData->Add(vmd);
			m_var_mesh_map[a_varLabel] = mesh3d;
			m_var_dim[a_varName] = 3;

			// also add bottom, surface and depth average state option
			avtVectorMetaData *vmd1 = new avtVectorMetaData(a_varLabel + m_surface_state_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd1);
			//debug1<<"add  "<<label+m_surface_state_suffix<<" ";
			avtVectorMetaData *vmd2 = new avtVectorMetaData(a_varLabel + m_bottom_state_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd2);
			//debug1<<"add  "<<label+m_bottom_state_suffix<<" ";
			avtVectorMetaData *vmd3 = new avtVectorMetaData(a_varLabel + m_depth_average_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd3);
			// debug1<<"add  "<< label+m_depth_average_suffix<<" ";
			m_var_name_label_map[a_varLabel + m_surface_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_bottom_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_depth_average_suffix] = a_varName;
			m_var_mesh_map[a_varLabel + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_depth_average_suffix] = m_mesh_2d;
		}
		else //3d scalar
		{
			avtScalarMetaData *smd = new avtScalarMetaData(a_varLabel, mesh3d, faceCent);
			a_metaData->Add(smd);
			m_var_mesh_map[a_varLabel] = mesh3d;

			// also add bottom, surface and depth average state option
			a_metaData->Add(new avtScalarMetaData(a_varLabel + m_surface_state_suffix, m_mesh_2d, faceCent));
			a_metaData->Add(new avtScalarMetaData(a_varLabel + m_bottom_state_suffix, m_mesh_2d, faceCent));
			a_metaData->Add(new avtScalarMetaData(a_varLabel + m_depth_average_suffix, m_mesh_2d, faceCent));
			//debug1<<"add  "<< label+m_depth_average_suffix<<" ";
			// all those surface, botootm and average are based on original data set
			m_var_name_label_map[a_varLabel + m_surface_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_bottom_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_depth_average_suffix] = a_varName;

			m_var_mesh_map[a_varLabel + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_depth_average_suffix] = m_mesh_2d;
			m_var_dim[a_varName] = 3;

		}
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, faceCent);
	}
	else
	{

	}


}
void    avtMDSCHISMFileFormat::addNodeCenterData(avtDatabaseMetaData * a_metaData,
	                                             SCHISMVar10            * a_varPtr,
	                                             const std::string   & a_varName,
	                                             const std::string   & a_varLabel,
	                                             const avtCentering  & a_center)
{

	std::string varName(a_varName);
	std::string label(a_varLabel);
	avtCentering avtCenter(a_center);

	//  scalar var on 2D
	if (a_varPtr->num_dims() <= 2)
	{
		a_metaData->Add(new avtScalarMetaData(label, m_mesh_2d, avtCenter));
		m_var_mesh_map[label] = m_mesh_2d;
		m_var_dim[varName] = a_varPtr->num_dims();
		debug1 << "added 2d scalar:" << label;
	}
	//  vector var on 2D
	else if ((a_varPtr->num_dims() == 3) && (!m_impl->SCHISMVarIs3D(a_varPtr)))
	{
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);
		a_metaData->Add(new avtVectorMetaData(label, m_mesh_2d, avtCenter, ucomps));
		m_var_mesh_map[label] = m_mesh_2d;
		m_var_dim[varName] = 2;

	}
	//  scalar var having layer dim
	else if ((a_varPtr->num_dims() == 3) && (m_impl->SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_mesh->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}

		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, m_mesh_3d, avtCenter);
		a_metaData->Add(new avtScalarMetaData(label, m_mesh_3d, avtCenter));
		m_var_mesh_map[label] = m_mesh_3d;
		// also add bottom, surface and depth average state option
		if (!(varName == MeshConstants10::ZCOORD))
		{
			a_metaData->Add(new avtScalarMetaData(
				label + m_surface_state_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_surface_state_suffix << " ";
			a_metaData->Add(new avtScalarMetaData(
				label + m_bottom_state_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_bottom_state_suffix << " ";
			a_metaData->Add(new avtScalarMetaData(
				label + m_depth_average_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_depth_average_suffix << " ";


			// all those surface, botootm and average are based on original data set
			m_var_name_label_map[label + m_surface_state_suffix] = varName;
			m_var_name_label_map[label + m_bottom_state_suffix] = varName;
			m_var_name_label_map[label + m_depth_average_suffix] = varName;

			m_var_mesh_map[label + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_depth_average_suffix] = m_mesh_2d;
		}
		m_var_dim[varName] = 3;

	}
	else if ((a_varPtr->num_dims() == 4) && (m_impl->SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_mesh->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		// last dim is vector component         
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);

		a_metaData->Add(new avtVectorMetaData(label, m_mesh_3d, avtCenter, ucomps));
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, m_mesh_3d, avtCenter);
		m_var_mesh_map[label] = m_mesh_3d;

		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtVectorMetaData(
			label + m_surface_state_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_bottom_state_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_depth_average_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;
		m_var_mesh_map[label + m_surface_state_suffix] = m_mesh_2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = m_mesh_2d;
		m_var_mesh_map[label + m_depth_average_suffix] = m_mesh_2d;
		m_var_dim[varName] = 3;

	}
	else
	{

	}
}

void    avtMDSCHISMFileFormat::addSideCenterData(avtDatabaseMetaData * a_metaData,
	SCHISMVar10           * a_varPtr,
	const std::string   & a_varName,
	const std::string   & a_varLabel,
	const avtCentering  & a_center)
{
	// only add centered  var now 
	string mesh2d = m_side_center_point_2d_mesh;
	string mesh3d = m_side_center_point_3d_mesh;
	std::string level_center = a_varPtr->get_vertical_center();

	if (level_center == MeshConstants10::HALF_LAYER)
	{
		mesh3d = m_face_center_point_3d_mesh;
	}

	avtCentering  nodeCent(AVT_NODECENT);
	std::string varName(a_varName);
	std::string label(a_varLabel);
	avtCentering avtCenter(a_center);

	//  scalar var on 2D
	if (a_varPtr->num_dims() <= 2)
	{
		a_metaData->Add(new avtScalarMetaData(label, mesh2d, avtCenter));
		m_var_mesh_map[label] = mesh2d;
		m_var_dim[varName] = a_varPtr->num_dims();
		debug1 << "added 2d scalar:" << label;
	}
	//  vector var on 2D
	else if ((a_varPtr->num_dims() == 3) && (!m_impl->SCHISMVarIs3D(a_varPtr)))
	{
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);
		a_metaData->Add(new avtVectorMetaData(label, mesh2d, avtCenter, ucomps));
		m_var_mesh_map[label] = mesh2d;
		m_var_dim[varName] = 2;

	}
	//  scalar var having layer dim
	else if ((a_varPtr->num_dims() == 3) && (m_impl->SCHISMVarIs3D(a_varPtr)))
	{

		if (!(m_mesh->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, avtCenter);
		a_metaData->Add(new avtScalarMetaData(label, mesh3d, avtCenter));
		m_var_mesh_map[label] = mesh3d;
		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtScalarMetaData(
			label + m_surface_state_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtScalarMetaData(
			label + m_bottom_state_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtScalarMetaData(
			label + m_depth_average_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		// all those surface, botootm and average are based on original data set
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;

		m_var_mesh_map[label + m_surface_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_depth_average_suffix] = mesh2d;
		m_var_dim[varName] = 3;

	}
	else if ((a_varPtr->num_dims() == 4) && (m_impl->SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_mesh->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		// last dim is vector component         
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);

		a_metaData->Add(new avtVectorMetaData(label, mesh3d, avtCenter, ucomps));
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, avtCenter);
		m_var_mesh_map[label] = mesh3d;

		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtVectorMetaData(
			label + m_surface_state_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_bottom_state_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_depth_average_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;
		m_var_mesh_map[label + m_surface_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_depth_average_suffix] = mesh2d;
		m_var_dim[varName] = 3;
	}
	else
	{
	}

}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

vtkDataSet *
avtMDSCHISMFileFormat::GetMesh(int timestate, int domainid, const char *meshname)
{

#ifdef PARALLEL
    int myrank = PAR_Rank();
	debug1 << "rank " << myrank << " getting mesh " << domainid << "\n";
#endif
    return m_impl->GetMesh(timestate,domainid,this,meshname);
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

vtkDataArray *
avtMDSCHISMFileFormat::GetVar(int timestate, int domain, const char *varname)
{

   return m_impl->GetVar(timestate,domain,varname);

   
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Fri Jul 31 10:54:14 PDT 2020
//
// ****************************************************************************

vtkDataArray *
avtMDSCHISMFileFormat::GetVectorVar(int timestate, int domain,const char *varname)
{

    return m_impl->GetVectorVar(timestate,domain,varname);
}


//avtVariableCache * avtMDSCHISMFileFormat::get_cache()
//{
//	return this->cache;
//}


